Collisions of charged black holes 
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We perform fully non-linear numerical simulations of charged-black-hole collisions, described by 
the Einstein-Maxwell equations, and contrast the results against analytic expectations. We focus 
on head-on collisions of non-spinning black holes, starting from rest and with the same charge 
to mass ratio, Q/M. The addition of charge to black holes introduces a new interesting channel 
of radiation and dynamics, most of which seem to be captured by Newtonian dynamics and flat- 
space intuition. The waveforms can be qualitatively described in terms of three stages; (i) an infall 
phase prior to the formation of a common apparent horizon; (ii) a nonlinear merger phase which 
corresponds to a peak in gravitational and electromagnetic energy; (iii) the ringdown marked by 
an oscillatory pattern with exponentially decaying amplitude and characteristic frequencies that are 
in good agreement with perturbative predictions. We observe that the amount of gravitational- 
wave energy generated throughout the collision decreases by about three orders of magnitude as the 
charge-to-mass ratio Q/M is increased from to 0.98. We interpret this decrease as a consequence 
of the smaller accelerations present for larger values of the charge. In contrast, the ratio of energy 
carried by electromagnetic to gravitational radiation increases, reaching about 22% for the maximum 
Q/M ratio explored, which is in good agreement with analytic predictions. 



I. INTRODUCTION 

In recent years, numerical relativity (NR) has gener- 
ated a wealth of information about astrophysical black- 
hole-binary systems; see [H-Q for the first complete sim- 
ulations and e.g. P4lCj| for a representative list of more 
recent studies. Results about the dynamics of black holes 
thus obtained are now actively employed in techniques 
and searches for gravitational wave signals in present and 
future generation gravitational wave detectors [TlT[l5j . 

While black-hole binaries interacting with electromag- 
netic fields and plasmas have been the subject of re- 
cent numerical studies (e.g. [D, @), the dynamics of bi- 
nary systems of charged, i.e. Reissner- Nordstrom (RN), 
black holes remain unexplored territory. Perhaps, this 
is due to the expectation that astrophysical black holes 
carry zero or very small charge; in particular, black 
holes with mass M, charge Q and angular momentum 
aM 2 are expected to discharge very quickly if Q/M > 
10- 13 (a/Af)- 1 / 2 (A//M Q ) 1 / 2 01 03. 

In spite of this expectation, however, there is a good 
deal of motivation for detailed investigations of the dy- 
namics of charged black holes. 

We first note that RN black holes possess a unique 
property amongst the black hole solutions of Einstein- 
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Maxwell theory in four dimensions. They possess an 
extremal limit which can be used to construct a static, 
regular (on and outside the event horizon) multi-black 
hole configuration fl8l (described by the Majumdar- 
Papapetrou solution |l9l. Ho|). This configuration can 
be interpreted as an exact cancellation, at each point, of 
the attractive (gravitational) and repulsive (electromag- 
netic) interactions — a no force condition. This condition 
is typically (but not always) associated to supersymmet- 
ric configurations and indeed the extremal RN solution is 
the only black hole solution in four dimensional Einstein- 
Maxwell theory that admits Killing spinors, when the 
theory is regarded as the bosonic sector of Af = 2 Super- 
gravity [2lL |22| . A natural question concerning the mod- 
elling of RN black holes in NR is how close can we get 
to extremality and hence consider the dynamics of these 
very special black holes. The ability to model such sys- 
tems could provide interesting applications. For instance 
it is possible to study analytically the dynamics of a per- 
turbed Majumdar-Papapetrou solution in the so-called 
moduli space approximation (23l . l2~i| . It would be inter- 
esting to compare this analytic approximation method 
with a fully non-linear NR simulation. 

Motivation for the numerical modelling of charged 
black holes also arises in the context of high energy col- 
lisions. It is expected that trans-Planckian particle col- 
lisions form black holes; moreover, well above the fun- 
damental Planck scale such processes should be well 
described by general relativity and other interactions 
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should become negligible [2 51 ], an idea poetically stated 
as matter does not matter for ultra high energy colli- 
sions . But is this expectation really correct? Calcu- 
lations of shock wave collisions suggest that even though 
other interactions — say charge — may become irrelevant 
in the ultra- relativistic limit, the properties of the final 
black hole (and of the associated emission of gravitational 
radiation) do depend on the amount of charge carried by 
the colliding particles [13, ■ This issue can be clari- 
fied by the simulation of high-energy collisions of charged 
black holes in the framework of NR and the subsequent 
comparison of the results to those obtained for electri- 
cally neutral systems. Recent works in this direction in- 
clude (29l - [35l | for binary black holes and [26| for boson 
stars. These, together with related incipient efforts to 
study gravity in higher-dimensional space-times [36l - l40j 
illustrate recent applications of numerical simulations to 
shed light on problems beyond astrophysical settings. 

In the context of astrophysics, charged black holes may 
be of interest in realistic systems. First, a rotating black 
hole in an external magnetic field will accrete charged 
particles up to a given value, Q = 2B J HH- Thus it 
is conceivable that astrophysical black holes could have 
some (albeit rather small) amount of electrical charge. 
Then it is of interest to understand the role of this charge 
in the Blandford-Znajck mechanism fl7j . which has been 
suggested for extracting spin energy from the hole, or in 
a related mechanism cap able of extracting energy from a 
moving black hole @,|42| to power outflows from accretion 
disk-fed black holes. NR simulations of charged black 
holes interacting with matter and surrounding plasma 
will enable us to study such effects. 

Finally we note a variety of conceptual aspects that 
merit a more detailed investigation of charged black-hole 
systems. In head-on collisions with small velocity, the 
intuition borrowed from Larmor's formula in Minkowski 
space suggests a steady growth of the emitted power with 
the acceleration. However, it is by now well established 
that for uncharged black holes the gravitational radiation 
strongly peaks near the time of formation of a common 
apparent horizon. Does the electromagnetic radiation 
emission follow a similar pattern? And what is the rel- 
ative fraction of electromagnetic to gravitational wave 
emissions? Moreover, a non-head on collision of charged 
non-spinning black holes will allow us to study, as the 
end state, a (perturbed) Kerr-Newman geometry, which 
would be extremely interesting: linearized perturbations 
around Kerr-Newman black holes do not decouple f43ll44j 
and so far close to nothing is known about their proper- 
ties. Among others, the stability of the Kerr-Newman 
metric is an outstanding open issue. Furthermore, it 
has been observed that the inspiral phase of an orbit- 
ing black-hole-binary system can be well understood via 
post- Newtonian methods [45| (see also e.g. (4(|[47j]). The 
additional radiative channel opened by the presence of 
electric charge provides additional scope to probe this 
observation. 

With the above motivation in mind we here initiate 



the numerical study of non- linear dynamics of binary sys- 
tems of RN black holes, building on previous numerical 
evolutions of the Einstein-Maxwell system [HI. l48l - [50| . For 
reasons of simplicity, we focus in this study on binary sys- 
tems for which initial data can be constructed by purely 
analytic means [5ll. l52l] : head-on collisions, starting from 
rest, of non-spinning black holes with equal chargc-to- 
mass ratio. This implies in particular that the black holes 
carry a charge of the same sign so that the electromag- 
netic force will always be repulsive. We will extract both 
gravitational and electromagnetic radiation and monitor 
their behaviour as the charge-to-mass-ratio parameter of 
the system is varied. 

For this purpose, we present in Sec. |TT] the evolution 
equations and the initial data used. In Sec. IIHI the 
method for extraction of gravitational and electromag- 
netic radiation is discussed. In Sec. IIVI we summarize 
our analytic calculations and compare in Sec. |V] their 
predictions with the numerical results. Throughout this 
work greek "spacetimc indices" run from to 3 and latin 
"spatial" indices from 1 to 3. 



II. EVOLUTION EQUATIONS 

In this paper we adopt the approach outlined in [49l.[53j 
to evolve the electro- vacuum Einstein-Maxwell equations 
which incorporates suitably added additional fields to en- 
sure the evolution will preserve the constraints. This 
amounts to considering an enlarged system of the form 



(2.1) 



where -kF^ v denotes the Hodge dual of the Maxwell- 
Faraday tensor F^ v , k is a constant and the four- 
velocity of the Eulerian observer. We recover the stan- 
dard Einstein-Maxwell system of equations when ^ = 
= $. With the scalar field 'J and pseudo-scalar $ in- 
troduced in this way, the natural evolution of this system 
drives ^ and $ to zero (for positive k), thus ensuring the 
magnetic and electric constraints are controlled [H, [53| . 
The electromagnetic stress-energy tensor takes the usual 
form 



T =J- 

47T 



Fp^Fux — —g^F^ F\ a 



(2.2) 



A. 3 + 1 decomposition 

We employ a Cauchy approach so we introduce a 3 + 1 
decomposition of all dynamical quantities. Concretely, 
we introduce the 3-mctric 
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and decompose the Maxwell-Faraday tensor into the 
more familiar electric and magnetic fields measured by 
the Eulcrian observer moving with four velocity n M 



F„ 



n v E^ + e^ va pB a n^ , 



(2.4) 



where we use the convention 61230 = \f~9, £a/3~f = 

tap-fgn 5 , £123 = y/j- 

Writing the evolution equations in the BSSN form (see, 
e.g., [Ill, |j5f| for details), we have, for the "gravitational" 
part 
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XHj 



X = 7 -V3 
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(dt-Cp)iij 

(dt-Cp)x 
(d t -Cp)K 

(d t -Cp)A ij = [■••]- 8ira ( X Sij " 

(d t -Cp)r = [■■■] -IQirax^f 



-2aAij , 

+ 47tq( ( o + 5) , 



- jk 5 

(2.6) 



where [• • • ] denotes the standard right-hand side of the 
BSSN equations in the absence of source terms. For the 
case of the electromagnetic energy-momentum tensor of 
Eqs. (|2.2[) , (|2.ip . the source terms are given by 
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p = T^n^ = — (E z 
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-EiEj - BiBj + \p l3 (E 2 + B 2 
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and S = j t3 Sij. The evolution of the electromagnetic 
fields is determined by Eq. (|2.ip whose 3+1 decomposi- 
tion becomes [501 



(dt - Cp) E i = aKE 1 + e ijk X ~ l [luB l d ]a + a (B l d^ kl + %id 3 B l - x'^klB' l d jX )] - a X f j d^ , 
(d t - Cp) B* = aKB* - e^x- 1 [%iE l d ja + a (E l drf m + judjE 1 - x^lkiE l d ]X )] - a X f'd^ , 



(d t - Cp) * = -aViE 1 - arf 



(d t -Cp)<f> = -aV t B l - a K $ 



Here, Cp denotes the Lie derivative along the shift vector 
/3\ The Hamiltonian and momentum constraint are 

n = 3 R + K 2 - K ij K l3 - 167rp = , 

M l = DjA l J - ^ASx^djX - \diK - Swji = , 

. (2.9) 

where Di is the covariant derivative associated with the 
three-metric 7^. 



B. Initial data 

As already mentioned in the Introduction, we focus 
here on black-hole binaries with equal charge and mass 
colliding from rest. For these configurations, it is pos- 
sible to construct initial data using the Brill-Lindquist 
construction [5l| (see also [13 )• The main ingredients of 
this procedure are as follows. 

For a vanishing shift f3' , time symmetry implies Kij = 
0. Combined with the condition of an initially vanishing 



magnetic field, the magnetic constraint DiB 1 = and 
momentum constraint are automatically satisfied. By 
further assuming the spatial metric to be conformally 
flat 

-,,d.r'd.r' = ^ (dx 2 + dy 2 + dz 2 ) , (2.10) 
the Hamiltonian constraint reduces to 

= -\e 2 iP 5 , (2.11) 

where A is the flat space Laplace operator. The electric 
constraint, Gauss's law, has the usual form 

DiE^Q. (2.12) 

Quite remarkably, for systems of black holes with equal 
charge-to-mass ratio, these equations have known ana- 
lytical solutions [52j ■ For the special case of two black 
holes momentarily at rest with "bare masses" mi, m-i 
and "bare charges" q±, (72 = 91^2/^1 this analytic solu- 
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tion is given by 



to the unit timclikc vector n M : 



mi 



m 2 



2 a? — xi\ 2\x — x-2\ 



1 ( qi 



<l2 



4 \\X — Xi\ \x — X2\ 



(2.13) 



X - .Tl| 



x - x 2 \ 



where Xi is the coordinate location of the ith "punc- 
ture" 

The initial data arc thus completely specified in terms 
of the independent mass and charge parameters mi, 1112, 
qi and the initial coordinate separation d of the holes. 
These uniquely determine the remaining charge param- 
eter q 2 via the condition of equal charge-to-mass ratio. 
In this study we always choose mi = m-i and, without 
loss of generality, position the two holes symmetrically 
around the origin such that Z\ = d/2 = —z 2 . The result- 
ing initial three metric 7^ follows from Eqs. (|2.10[) . (|2.13[) 
while the extrinsic curvature Kij and magnetic field B l 
vanish on the initial slice. 

Finally, the time evolution of the fields is determined 
by Eqs. (|2.6p and (|2.8[) . We use the same gauge con- 
ditions and outer boundary conditions for the BSSN 
variables as used in vacuum simulations [5fj|. As outer 
boundary condition for the electric and magnetic fields 
we have imposed a falloff as 1/r 2 — from (|2.13[) . For the 
additional scalar fields a satisfactory behaviour is ob- 
served by imposing a falloff as 1 /r 3 (which is the expected 
falloff rate from dimensional grounds). 



III. WAVE EXTRACTION 

For a given set of initial parameters mi = m2, qi =92, 
d, the time evolution provides us with the spatial met- 
ric "fij , the extrinsic curvature as well as the electric 
and magnetic fields E % , B % as functions of time. These 
fields enable us to extract the gravitational and electro- 
magnetic radiation as follows. 

For the gravitational wave signal we calculate the 
Newman-Penrose scalar ^4 defined as 



*4 = C a p lS k a fh f) k' y m s 



(3.1) 



where C a p.ys is the Weyl tensor and k 7 m are part of a 
null tetrad l,k 7 m,m satisfying — I ■ k = 1 = m ■ m; all 
other inner products vanish. In practice I, k and m are 
constructed from an orthonormal triad u, V, W orthogonal 



1 We note that this foliation, in isotropic coordinates, only covers 
the outside of the external horizon. 



r 
k° 



-j={n a +u a ) , 

-L(„« + ,w<) . 



(3.2) 



We refer the interested reader to |57j for more details 
about the numerical implementation and [58j for a re- 
view of the formalism; here we merely note that asymp- 
totically the triad vectors u, v, w behave as the unit 
radial, polar and azimuthal vectors r, 8, (f>. 

Similarly, we extract the electromagnetic wave signal in 
the form of the scalar functions, $1 and $2 [HI, defined 
as 

1 



$1 = -F^ {n u + m^m") 



(3.3) 

$ 2 = F^mW . (3.4) 
For outgoing waves at infinity, these quantities behave as 



- (E f + iB f ) 



Eg - iE x 



(3.5) 



At a given extraction radius i?ox, we perform a multi- 
polar decomposition by projecting ^4, <!>i and $2 onto 
spherical harmonics of spin weight s = —2, and — 1 
respectively: 



I . n 1 



<J>l(M^) = E^™^ 



I . n 1 



$ 2 (t,9,<t>) = J24 m W l - m 1 (9,<l>) 

I . rj 1 



(3.6) 
(3.7) 
(3.8) 



In terms of these multipoles. the radiated flux and energy 
is given by the expressions |59j 



Fqw — 



dE, 



GW 



dt 



lim T^E f dt'i> lm (t') 



EM 



dEvM 
dt 



r— >oc 47T — * 



(3.9) 
(3.10) 



Lm 



As is well known from simulations of uncharged black- 
hole binaries, initial data obtained from the Brill- 
Lindquist construction contain "spurious" radiation, 
which is an artifact of the conformal-flatncss assumption. 
In calculating properties of the radiation, we account for 
this effect by starting the integration of the radiated flux 
in Eqs. p.9[) , (|3 . 10[) at some finite time At after the start 
of the simulation, thus allowing the spurious pulse to first 
radiate off the computational domain. In practice, we ob- 
tain satisfactory results by choosing At = i? ex + 50 M. 
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Because the physical radiation is very weak for both the 
gravitational and electromagnetic channel in this early 
infall stage, the error incurred by this truncation is neg- 
ligible compared with the uncertainties due to discretiza- 
tion; cf. Sec. IVDl 



IV. ANALYTIC PREDICTIONS 

Before discussing in detail the results of our numerical 
simulations, it is instructive to discuss the behaviour of 
the binary system as expected from an analytic approx- 
imation. Such an analysis not only serves an intuitive 
understanding of the binary's dynamics, but also pro- 
vides predictions to compare with the numerical results 
presented below. 

For this purpose we consider the electrodynamics of a 
system of two equal point charges in a Minkowski back- 
ground spacetime. As in the black-hole case, we denote 
by 9i = 92 = Q/2 and mi = to 2 = M/2 the electric 
charge and mass of the particles which are initially at 
rest at position z = ±d/2. 

It turns out to be useful to first consider point charges 
in Minkowski spacetime in the static limit. The expected 
behaviour of the radial component of the resulting elec- 
tric field is given by [6(| 



^ l + l Y lm {e,y) 
Ef = 4tt y j nl t A qi r 



I, lit 



21 + 1 



(4.1) 



which for a system of two charges of equal magnitude at 
z = ±d/2 becomes 



Ef ~ VAttQ— + \ —Qd —r 
r z V 20 r 4 



(4.2) 



The dipole vanishes in this case due to the reflection 
symmetry across z = 0. This symmetry is naturally 
preserved during the time evolution of the two-charge 
system. Furthermore, the total electric charge Q is con- 
served so that the leading-order behaviour of the electro- 
magnetic radiation is given by variation of the electric 
quadrupole, just as for the gravitational radiation. No- 
tice that in principle other radiative contributions can 
arise from the accelerated motion of the charged black 
holes. From experience with gravitational radiation gen- 
erated in the collision of electrically neutral black-hole 
binaries, however, we expect this "Bremsstrahlung" to be 
small in comparison with the merger signal and hence ig- 
nore its contributions in this simple approximation. The 
good agreement with the numerical results presented in 
the next section bears out the validity of this quadrupole 
approximation. In consequence, it appears legitimate to 
regard the "strength" of the collision and the excitation 
of the black-hole ringdown to be purely kinematic effects. 

An estimate for the monopole and quadrupole ampli- 
tudes in the limit of two static point charges is then ob- 
tained from inserting the radial component of the electric 



field (14.21) into the expression (|3.5|) for $1 and its multi- 
polar decomposition (|3.7[) 

r^f = v^Q 



r^f 



1.77Q, 

0.59Qd 2 



(4.3) 
(4.4) 



The expectation is that these expressions provide a good 
approximation for the wave signal during the early infall 
stage when the black holes are moving with small veloc- 
ities. Equation (|4.3|) should also provide a good approx- 
imation for (j>i° after the merger and ringdown whereas 
the quadrupole cj) 20 should eventually approach zero as 
a single merged hole corresponds to the case d = in 
Eq. dill) . 

To obtain analytic estimates for the collision time and 
the emitted radiation, we need to describe the dynamic 
behaviour of the two point charges. Our starting point 
for this discussion is the combined gravitational and elec- 
tromagnetic potential energy for two charges i = 1, 2 in 
Minkowski spacetime with mass and charge m, , % at dis- 
tance r from each other 



V=- 



1 



<7i<72 



47re r 



(4.5) 



For the case of two charges with equal mass and charge 
rrii = M/2, qi = Q/2 and starting from rest at z = 
id/2, conservation of energy implies 



Mz 2 



MB 



MB 



(4.6) 



4z 2d ' 

where we have used units with G = 47reo = 1 and 

B=l- Q 2 /M 2 . (4.7) 

The resulting equation of motion for z(t) is obtained by 
differentiating Eq. (|4.6[) which results in 



Mz 



8z 2 
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= —M 



B 



(4.8) 



An estimate for the time for collision follows from inte- 
grating Eq. (|4~6|) over z e [d/2,0] 



(^collision \ 7T d 

M J ~ 2 3 M 3 B 



(4.9) 



From the dynamic evolution of the system we can 
derive an approximate prediction for the electromag- 
netic radiation by evaluating the (traceless) electric 
quadrupole tensor Qy = J d 3 xp(x)(3xiXj — r 2 Sij) [60l |. 
In terms of this quadrupole tensor, the total power radi- 
ated is given by [60| 

For clarity we have reinstated the factors 4ire and c 5 
here. Using 



dt 3 



(z 2 ) = 6zz + 2z'z , 



(4.11) 
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and the equations of motion (I4.6|) . (|4.8[) we find Using J dt(. . .) = J dz/z(. . .), we can evaluate the time 

integral up to some cutoff separation, say z m ; n = a b 6, 
where b is the horizon radius of the initial black hole, 

B 3 M 3 2 (l/ — 2/d) ^ = v ^)/^ ana - ah = ^(l) i s a constant. This 

F em = ^A? LL . (4.12) gives, 



1920z 4 

I 



^fad _ R 5/2 . ,3/2 p 2 (j ~ ^ b b)^ (15d 2 + 24ri Qb b + 32^g6 2 ) 

M _ ^ Q 50400(da 6 b)V2 ■ (4 - idj 

Emission of gravitational radiation follows from the quadrupole formula, which is a numerical factor 4 times larger, 
and where the charge is be replaced by the mass, 

ggZ _ .5/2 , f7 /2 (j ' 2a b b)3/^(15rf 2 + 24d« b fr + 32cgb 2 ) 
M 12600(da fe 6)7/2 ' ^ a4J 



For Q = 0, at = 1, d = oo we thus obtain 



M 



83o ^ °- 0012 ■ 



(4.15) 



in agreement to within a factor of 2 with numerical sim- 
ulations (see [33| and Table U below; the agreement could 
be improved by assuming a b ~ 1.3). As a general result 
of this analysis we find in this approximation, 



Q 2 

E GW 4M 2 



(4.16) 



For non-extremal holes Q < M, our analytic considera- 
tions therefore predict that the energy emitted in electro- 
magnetic radiation is at most 25% of the energy lost in 
gravitational radiation. As we shall see below, this turns 
out to be a remarkably good prediction for the results 
obtained from fully numerical simulations. 



V. NUMERICAL RESULTS 

The numerical integration of the Einstein-Maxwell 
equations (|2.6[) . (|2.8[) has been performed using fourth- 
order spatial discretization with the Lean code, origi- 
nally presented in (5?J for vacuum spacetimes. Lean 
is based on the Cactus Computational toolkit foil ], 
the Carpet mesh refinement package [62|, [63[ and uses 
AHFinderDirect for tracking apparent horizons [(34|, 
[65| . For further details of the numerical methods see 
Ref. [13. 

The initial parameters as well as the grid setup and the 
radiated gravitational and electromagnetic wave energy 
for our set of binary configurations is listed in Table HI 
All binaries start from rest with a coordinate distance 
d/M ~ 8 or d/M ~ 16 while the charge-to-mass ratio 
has been varied from Q/M = to Q/M = 0.98. Note 
that identical coordinate separations of the punctures for 
different values of the charge Q/M correspond to dif- 
ferent horizon-to-horizon proper distances. This differ- 
ence is expected and in fact analysis of the RN solution 



predicts a divergence of the proper distance in the limit 
Q/M -> 1. 



A. Code tests 

Before discussing the obtained results in more detail, 
we present two tests to validate the performance of our 
numerical implementation of the evolution equations, (i) 
Single black-hole evolutions in geodesic slicing which is 
known to result in numerical instabilities after relatively 
short times but facilitates direct comparison with a semi- 
analytic solution and (ii) Convergence analysis of the ra- 
diated quadrupole waveforms for simulation d08q05 of 
Table H 

The geodesic slicing condition is enforced by setting 
the gauge functions to a = 1, ft 1 = throughout the 
evolution. The space part of the Reissner-Nordstrom so- 
lution in isotropic coordinates is given by Eq. (|2 . 10[) with 
a conformal factor [69, k37| 



M_y qi_ 



(5.1) 



The time evolution of this solution is not known in closed 
analytic form, but the resulting metric components can 
be constructed straightforwardly via a simple integration 
procedure, cf. Appendix [S] As expected, we find a time 
evolution in this gauge to become numerically unstable 
at times r of a few M. Before the breaking down of the 
evolution, however, we can safely compare the numerical 
and "analytical" solutions. This comparison is shown in 
Fig. [T] for the j zz component of the spatial metric and 
the E z component of the electric field and demonstrates 
excellent agreement between the semi-analytic and nu- 
merical results. 

For the second test, we have evolved model d08q05 us- 
ing three different resolutions as listed in Table U and ex- 
tracted the gravitational and electromagnetic quadrupole 
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TABLE I. Grid structure in the notation of Sec. II E of [57}], coordinate distance d/M, proper horizon-to-horizon distance L/M, 
charge Q/M, gravitational (E^^) and electromagnetic (E™f) radiated energy for our set of simulations. The radiated energy 
has been computed using only the I = 2, m = mode; the energy contained higher-order multipoles such as / = 4, m — is 
negligible for all configurations. 



Run 


Grid 


d/M 


L/M 


Q/M 


TTi GW 




77EM 




d08q00 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


11.56 





5.1 X 10" 


4 


_ 




d08q03 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


11.60 


0.3 


4.5 x 10" 


-4 


1.3 x 10" 


5 


d08q04 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


11.65 


0.4 


4.0 x 10" 


-4 


2.1 x 10" 


-5 


d08q05c 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/64} 


8.002 


11.67 


0.5 


3.3 x 10" 


-4 


2.7 x 10" 


5 


d08q05m 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


11.70 


0.5 


3.4 x 10" 


-4 


2.7 x 10" 


5 


d08q05f 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/96} 


8.002 


11.67 


0.5 


3.4 x 10" 


-4 


2.7 x 10" 


-5 


d08q055 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


11.70 


0.55 


3.0 x 10" 


-4 


2.89 x 10 


-S 


d08q06 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


11.75 


0.6 


2.6 x 10" 


-4 


2.97 x 10 


-5 


d08q07 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


11.87 


0.7 


1.8 x 10" 


-4 


2.7 x 10" 


-5 


d08q08 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


12.0 


0.8 


9.8 x 10" 


-5 


1.8 x 10" 


-5 


d08q09 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


12.3 


0.9 


2.6 x 10" 


-5 


5.5 x 10" 


6 


d08q098cc 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/64} 


8.002 


12.3 


0.98 


7.0 x 10" 


7 


2.1 x 10" 


7 


d08q098c 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/80} 


8.002 


13.1 


0.98 


4.3 x 10" 


7 


1.4 x 10" 


7 


d08q098m 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/96} 


8.002 


13.1 


0.98 


3.4 x 10" 


7 


1.0 x 10" 


7 


d08q098f 


{(256,128,64,32,16,8) x (2, 1, 0.5), 1/112} 


8.002 


13.0 


0.98 


4.0 x 10" 


7 


9.5 x 10" 


-8 


d08q098ff 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/128} 


8.002 


13.0 


0.98 


4.05 x 10 


-7 


8.75 x 10 


-8 


d08q098fff 


{(256, 128, 64, 32, 16, 8) x (2, 1, 0.5), 1/136} 


8.002 


13.1 


0.98 


3.73 x 10 


-7 


8.41 x 10 


-8 


dl6q00 


{(256, 128, 64, 32, 16) x (4, 2, 1, 0.5), 1/64} 


16.002 


20.2 





5.5 x 10" 


4 






dl6q05 


{(256, 128, 64, 32, 16) x (4, 2, 1, 0.5), 1/64} 


16.002 


20.3 


0.5 


3.6 x 10" 


4 


2.9 x 10" 


-5 


dl6q08 


{(256, 128, 64, 32, 16) x (4, 2, 1, 0.5), 1/80} 


16.002 


20.7 


0.8 


1.05 x 10 


-4 


1.9 x 10" 


-5 


d!6q09 


{(256, 128, 64, 32, 16) x (4, 2, 1, 0.5), 1/80} 


16.002 


21.0 


0.9 


2.7 x 10" 


-5 


5.9 x 10" 


-6 




(/ = 2,to = 0) at R cx = 100 M. For fourth-order con- 
vergence, we expect the differences between the higher 
resolution simulations to be a factor 2.78 smaller than 
their coarser resolution counterparts. The numerically 
obtained differences are displayed with the corresponding 
rescaling in Fig. [2] Throughout the physically relevant 
part of the waveform, we observe the expected fourth- 
order convergence. Only the spurious initial radiation 



(cf. the discussion at the end of Sec. IIII[) at early times 
At < —20 in the figure exhibits convergence closer to sec- 
ond order, presumably a consequence of high-frequency 
noise contained in this spurious part of the signal. From 
Richardson extrapolation of our results we estimate the 
truncation error of the radiated waves to be about 1 %. 
The error due to extraction at finite radius, on the other 
hand, is estimated to be 2 % at i? GX = 100 M. 



- - 20 40 60 80 100 120 20 40 60 80 100 120 

A A 

FIG. 2. Convergence analysis for simulation d08q05 of Table [J with resolutions h c = M/64, h m = M/80 and hf = Af/96. The 
panels show differences of the (2, 0) multipoles of the real parts of ^4 (left) and $2 (right) extracted at i? GX = 100 M; in each 
case, the high-resolution differences have been rescaled by a factor 2.78 as expected for fourth-order convergence. 



B. Collisions of two black holes: the "static" 
components and infall time 

We start the discussion of our results with the be- 
haviour of the gravitational and electromagnetic multi- 
poles when the system is in a nearly static configura- 
tion, i.e. shortly after the start of the simulation and 
at late stages after the ringdown of the post-merger 
hole. At these times, we expect our analytic predictions 
(|4.3j) . (|4.4|) for the monopole and dipole of the electro- 
magnetic field to provide a rather accurate description. 
Furthermore, the total spacetimc charge Q is conserved 
throughout the evolution, so that the monopole compo- 
nent of $1 should be described by (|4.3p at all times. The 
quadrupole, on the other hand, is expected to deviate sig- 
nificantly from the static prediction (|4.4[) when the black 
holes start moving fast. 

As demonstrated in Fig. [3J we find our results to be 
consistent with this picture. Here we plot the monopole 
and quadrupole of $1. The monopole part (left panel) 
captures the Coulomb field and can thus be compared 
with the total charge of the system. It is constant 
throughout the evolution to within numerical error and 
shows agreement with the analytic prediction of Eq. (|4.3[) 
within numerical uncertainties; we measure a slightly 
smaller value for the monopole field than expected from 
the total charge of the system, but the measured value 
should increase with extraction radii and agree with the 
total charge expectation at infinity. This is consistent 
with the extrapolation of the measured value to infin- 
ity as shown in the figure. The quadrupole part (right 
panel) starts at a non-zero value in excellent agreement 
with Eq. (|4.4[) . deviates substantially during the highly 
dynamic plunge and merger stage and eventually rings 
down towards the static limit = as expected for a 
spherically symmetric charge distribution. 

The analytic approximation of Sec. IIVI also predicts 



a value for the time of collision (|4.9|) for a given set of 
initial parameters. In particular, we see from this predic- 
tion that for fixed initial separation d and mass M the 
collision time scales with the charge as ^collision ~ 1/vB. 
In comparing these predictions with our numerical re- 
sults we face the difficulty of not having an unambigu- 
ous definition of the separation of the black holes in the 
fully general relativistic case. From the entries in Table U 
we see that the proper distance L varies only mildly for 
fixed coordinate distance d up to Q/M f=s 0.8. For nearly 
extremal values of Q, however, L starts increasing sig- 
nificantly as expected from our discussion at the start of 
this section. We therefore expect the collision time of the 
numerical simulations rescaled by ^/B/to, where to is the 
corresponding time for the uncharged case, to be close to 
unity over a wide range of Q/M and show some devia- 
tion close to Q/M = 1. This expectation is borne out in 
Fig. [4] where we show this rescaled collision time, deter- 
mined numerically as the first appearance of a common 
apparent horizon, as a function of Q/M . 



C. Waveforms: infall, merger and ringdown 

The dynamical behaviour of all our simulations is qual- 
itatively well represented by the waveforms shown in 
Fig. [5] for simulations dl6q00, dl6q05 and dl6q09. The 
panels show the real part of the gravitational (left) and 
electromagnetic (right) quadrupole extracted at R ex = 
100 M as a function of time with At = defined as the 
time of the global maximum of the waveform. From the 
classical analysis (|4.10p . we expect the waveforms ^4, $2 
to scale roughly with B and the mass or charge of the 
black holes (the scaling with B is non-trivial, but both 
an analytic estimate and the numerical results indicate 
the scaling is approximately linear, which we shall there- 
fore use for re-scaling the plots in the figure). 
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300 




50 100 150 200 250 50 100 150 200 250 300 



t/M t/M 

FIG. 3. Monopole (j>i° (left) and quadrupole cf>i° (right) of the radial part of the electromagnetic field $i extracted at 7? ex = 
100 M for simulation d08q05 of Table [U The dashed curves show the predictions of Eqs. Q4.3[) . (|4.4[) at R = oo in the static 
limit. For the monopole case, we also added the curves obtained by extrapolating the results to infinite extraction radius; these 
curves — dotted lines — essentially overlap with the predictions from Eq. (|4.3[) . 



1.03 
1.02 
1.01 



0.99 
0.98 



' 0.0 0.2 0.4 0.6 0.8 1.0 

Q/M 

FIG. 4. Time for apparent horizon formation, re-scaled by 
the factor \[~B and the apparent horizon formation time to for 
an electrically neutral binary. We note that the change in the 
quantity we plot is only, at most, of 2%. The coordinate time 
itself, however, varies by a factor 5 as one goes from Q = to 
Q = 0.98M. 



The early stages of the signals are marked by the spu- 
rious radiation due to the construction of initial data 
which we ignore in our analysis. Following a relatively 
weak phase of wave emission during the infall of the 
holes, the radiation increases strongly during the black- 
hole merger around At = in the figure and decays ex- 
ponentially as the final hole rings down into a stationary 
state. This overall structure of the signals is rather sim- 
ilar for the electromagnetic and the gravitational parts 
and follows the main pattern observed for gravitational- 
wave emission in head-on collisions of uncharged black 

holes mm. 

The final, exponentially damped ringdown phase is 



Q/M 


QNM 


ext 





0.374 - 0.0890i 


0.374 - 0.088i 




0.458 - 0.0950; 




0.3 


0.376 - 0.0892; 


0.375 - 0.092i 




0.470 - 0.0958; 


o.48i - o.ioo; 


0.5 


0.382 - 0.0896; 


0.381 - 0.091; 




0.494 - 0.0972; 


0.511 - 0.096; 


0.9 


0.382 - 0.0896; 


0.381 - 0.091; 




0.494 - 0.0972; 


? 



TABLE II. Comparison of the ringdown frequencies obtained 
from (i) perturbative calculations [44| and (ii) fitting a two- 
mode profile to the numerically extracted waveforms. For 
Q/M = the electromagnetic modes are not excited. For 
values of Q/M > 0.9 the electromagnetic mode becomes so 
weak that we can no longer unambiguously identify it in the 
numerical data. 

well described by perturbation techniques [44[ . In partic- 
ular, charged black holes are expected to oscillate with 
two different types of modes, one of gravitational and 
one of electromagnetic origin. For the case of vanish- 
ing charge, the electromagnetic modes are not present, 
but they generally couple for charged black holes, and 
we expect both modes to be present in the spectra of our 
gravitational and electromagnetic waveforms. For verifi- 
cation we have fitted the late-stages of the waveforms to 
a two-mode, exponentially damped sinusoid waveform 

f(t) = A ie - luJlt + A 2 e^ 2 \ (5.2) 

where Ai are real-valued amplitudes and complex fre- 
quencies. The results are summarized in Table |TT] for 
selected values of the chargc-to-mass ratio of the post- 
merger black hole. Real and imaginary parts of the fitted 
frequencies agree within a few percent or better with the 
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FIG. 5. Real part of the (2, 0) mode of * 4 (left) and $ 2 (right) extracted at i? ox = 100M. 



perturbative predictions. For the large value Q/M, how- 
ever, the wave signal is very weak and in such good agree- 
ment with a single ringdown mode (the gravitational one) 
that we cannot clearly identify a second, electromagnetic 
component. This feature is explained once we understand 
how the total radiated energy is distributed between the 
gravitational and the electromagnetic channels. For this 
purpose, we plot in Fig. [6] the Fourier spectrum of the 
relevant wavefunctions or, more precisely, their dominant 
quadrupolc contributions obtained for simulation d08q03 
|</> 20 | 2 , 1'*/' 20 ! 2 , where for any function / 

/oo 
e^f(t)dt. (5.3) 
-oo 

It is clear from the figure that most of the energy is 
carried in the fundamental gravitational-wave like mode 



with a peak at approximately uj ~ 0.37, close to the os- 
cillation frequency of the fundamental gravitational ring- 
down mode; see Table HT1 

D. Radiated energy and fluxes 

The electromagnetic and gravitational wave fluxes are 
given by Eqs. (|3.9|) and (|3.10|) . We have already noticed 
from the waveforms in Fig. [5] that the electromagnetic 
signal follows a pattern quite similar to the gravitational 
one. The same holds for the energy flux which is shown 
in Fig. [7] for a subset of our simulations with Q/M = 0, 
0.5 and 0.9. From the figure, as well as the numbers in 
Table U we observe that the energy carried by gravita- 
tional radiation decreases with increasing Q/M, as the 
acceleration becomes smaller and quadrupole emission is 




FIG. 6. Power spectrum for the gravitational (long dashed) 
and electromagnetic (short dashed) quadrupole extracted 
from simulation d08q03. Note that the spectrum peaks 
near the fundamental ringdown frequency of the gravitational 
mode; cf. Table HH 



FIG. 7. Radiated fluxes for simulations d08q05, d08q09 and 
d08q00 of Table [U We have aligned the curves in time such 
that their global maximum coincides with t = 0. The inset 
shows the exact same plot with the y-axis in logarithmic units. 
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0.6 



0.5 



0.4 



0.3 



0.2 



0.1 



Radiated energy 



7 GW 
J rad 



EM /ciGW 



0.2 



0.4 0.6 

Q/M 



0.8 



FIG. 8. Energy radiated in the gravitational and electromag- 
netic quadrupole as well as the ratio of the two as a function 
of Q/M. 



suppressed, in agreement with prediction (|4.14|) . 

This is further illustrated in Fig. [HJ which illus- 
trates the radiated energy carried in the gravitational 
quadrupole and the electromagnetic quadrupole as well 
as their ratio as functions of the chargc-to-mass ratio 
Q/M. For the case of vanishing charge, the total radi- 
ated energy is already known from the literature; e.g. [33jj. 
The value increases mildly with the initial separation as 
a consequence of the slightly larger collision velocity but 
is generally found to be close to £gj /M = 0.055 %. Our 
values of 0.051 % for d/M ~ 8 and 0.055 % for d/M ~ 16 
are in good agreement with the literature. As we in- 
crease Q/M, however, E^± decreases significantly and 
for Q/M = 0.9 (0.98) has dropped by a factor of about 
20 (10 3 ) relative to the uncharged case. For practical 
reasons, we have explored the largest ratio Q/M = 0.98 
for the smaller initial separation d/M ~ 8 only; the near 
cancellation of the gravitational and electromagnetic in- 
teraction and the resulting slow-down of the collision lead 
to a very long infall stage with essentially zero dynamics. 

In contrast to the monotonically decreasing 
gravitational-wave energy, the electromagnetic sig- 
nal reaches a local maximum around Q/M = 0.6, an 
expected observation as the electromagnetic radiation 
necessarily vanishes for Q/M = (no charge) and 
Q/M = 1 (no acceleration) but takes on non-zero 
values in the regime in between. Closer analysis of 
our classical, flat-space calculation (|4.13j) predicts a 
maximum electromagnetic radiation output at 



Qv 



13 



14 



M « 0.605AJ, 



(5.4) 



in excellent agreement with the results of our NR simu- 
lations. 

We finally consider the ratio of electromagnetic to 
gravitational wave energy (dotted curve in Fig. |5). As 
predicted by our analytic calculation (|4.16[) . this ratio 



increases monotonically with Q/M for fixed separation 
d. A fit of our numerical results yields E™/E^ = 
0.27 Q 2 /M 2 and for our largest value Q/M = 0.98, we 
obtain a ratio of 0.227 to be compared with ~ 0.24 as 
predicted by Eq. (|4. 16|) . Bearing in mind the simplicity 
of our analytic model in Sec. IIV1 the quantitative agree- 
ment is remarkable. 



VI. FINAL REMARKS 

We have performed a numerical study of collisions of 
charged black holes with equal mass and charge in the 
framework of the fully non- linear Einstein-Maxwell equa- 
tions. Our first observation is that the numerical relativ- 
ity techniques (formulation of the evolution equations, 
gauge conditions and initial data construction) devel- 
oped for electrically neutral black-hole binaries can be 
straightforwardly extended to successfully model charged 
binaries even for nearly extremal charge-to-mass ratios 
Q/M < 1. In particular, we notice the contrast with the 
case of rotating black holes with nearly extremal spin 
which represents a more delicate task for state-of-the-art 
numerical relativity; cf. Refs. [gl, [69| for the latest de- 
velopments on this front. This absence of difficulties for 
charged holes is not entirely unexpected. Considering the 
construction of initial data, for instance, an important 
difference arises in the customary choice of conformally 
flat Bowen-York initial data [7(| which greatly simplifies 
the initial data problem. While the Kerr solution for 
a single rotating black hole does not admit conformally 
flat slices [tT| and therefore inevitably results in spurious 
radiation, especially for large spin parameters, this diffi- 
culty does not arise for charged, but non-rotating black 
holes; cf. Eq. (JO) and [H]. 

The excellent agreement between the classical calcu- 
lation for the energy emission and the numerical results 
reported here, allow for an investigation of Cosmic Cen- 
sorship close to extremality. If we take two black holes 
with Mi = M 2 = M/2, Q 1 = Q 2 = (M - S)/2 and we 
let them fall from infinity, to first order in <5 we get 



Qtot = M — 8 
M,^ = M - E Ta . d 



(6.1) 



Now, the classical result (|4.14[) implies that the domi- 
nant term for the radiated energy is -E ra d ~ B 5 ^ 2 M ~ 
(S/M) 5 / 2 M . Thus we get 



Qi 



Aft, 



1 



S 

M 



5/2 



(6.2) 



where k is a constant. We get the striking conclusion that 
Cosmic Censorship is preserved for charged collisions of 
nearly extremal holes (6 <C Af), on account of the much 
longer collision time, which yields much lower velocities 
and therefore much lower energy output. The differences 
between the cases of spinning mergers and charged colli- 
sions are interesting. In the former case, naked singulari- 
ties are avoided by radiation carrying away more angular 
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momentum (via orbital hangup [72[). In the latter case, 
our results suggest that naked singularities are avoided 
by the smaller radiation emission, due to the smaller ac- 
celerations involved in the infall. 

It is even possible to construct binary initial data in 
closed analytic form, analogous to that of Brill-Lindquist 
data, for the special case of non-spinning binaries with 
equal charge-to-mass ratio starting from rest and we 
have restricted our present study to this case. Specifi- 
cally, we have evolved a sequence of binaries with Q /M 
varying from zero to values close to extremality. Start- 
ing with the electrically neutral case, where our gravita- 
tional wave emission E^W /M = 0.055 % agrees well with 
the literature, we observe a monotonic decrease of the 
emitted gravitational wave energy as we increase Q/M. 
For our largest value Q/M = 0.98, E^Jf is reduced by 
about three orders of magnitude, as the near cancella- 
tion of the gravitational and electromagnetic forces sub- 
stantially slows down the collision. In contrast, the ra- 
diated electromagnetic energy reaches a maximum near 
Q/M = 0.6 but always remains significantly below its 
gravitational counterpart. Indeed, the ratio Ef^/E^^ 
increases monotonically with Q/M and approaches about 
25 % in the limit Q/M -» 1. We find all these results 
to be in remarkably good qualitative and quantitative 
agreement with analytic approximations obtained in the 
framework of the dynamics of two point charges in a 
Minkowski background. This approximation also pre- 
dicts that the collision t ime relative to that of the un- 
charged case scales ~ y/l — Q 2 /M 2 which is confirmed 
within a few percent by our numerical simulations. 

Our present study paves the way for various future ex- 
tensions. Quite naturally, it will be important to consider 
more generic types of initial data in order to tackle some 
of the issues discussed in the Introduction. A non-zero 
boost, for instance, will allow us to study both binary 
black hole systems that will coalesce into a Kerr- Newman 
black hole and the impact of electric charge on the dy- 
namics of wave emission (electromagnetic and gravita- 
tional) in high energy collisions. In this context the ro- 
bustness of our simulations is particularly encouraging, 
as we have not encountered stability issues as observed 
in the study of black- hole collisions in higher-dimensional 
spacetimes [73| . 

A further interesting extension presently under study 
is the case of oppositely charged black holes. Quite likely, 
the remarkable accuracy of our simple analytic models is 
in part due to the relatively small, "non-rclativistic" col- 
lision speeds caused by the electric repulsion of the equal 
charges. Furthermore, the gravitational quadrupole for- 
mula (|4.14[) will still apply for opposite charges, but then 
B = 1 + Q 2 /M 2 , and the formula predicts an enhance- 
ment of almost two orders of magnitude in the gravi- 
tational radiation emitted when going from Q = to 
Q = M (without accounting for additional contributions 
due to dipole radiation and to "Brcmsstrahlung" by ac- 
celerated charges). This would release about 3% of the 
total center of mass energy as gravitational waves. Even 



more impressive is the possibility of observing a huge 
splash of electromagnetic energy when both holes are 
nearly extremal. The area theorem, which yields a poor 
estimate in the neutral case, bounds the total radiation 
to be less than ~ 65% the CM energy; how close one gets 
to this number is up to nonlinear evolutions of the kind 
reported in this work. 
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Appendix A: Geodesic slicing 

In the usual Schwarzschild-like coordinates, the 
Reissncr-Nordstrom line element and electromagnetic po- 
tential are given by 
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ds 2 



-f(R)dt 2 



W) 



(Al) 



where f(R) = 1 — ^ + fp. For a radially in- falling 
massive particle, starting from rest at R = Rq, the energy 
per unit mass is y/ f(R ). The geodesic equation (for in- 
falling particles) may then be written as 



f = W- 7F~>/K5Fm ( A2) 

With these equations and the initial condition R(t = 
0) = Rq, we can numerically integrate this system and 
thus have R = R(t,Rq). Assuming such a coordinate 



transformation, (t,R) (r, Rq), the metric takes the 
form 



ds" = —dr 



(A3) 

It remains now to perform the coordinate transformation 
Rq — > r that guarantees the metric an isotropic form at 
t = 0. This can be accomplished with 



dr 



r R ^j%Roj' 
integrating we obtain 

M 



Mr) 



1 



2 Q 2 



2r J Ar 2 
The final form for the metric is then 



(A4) 



(A5) 



ds 2 = -dr 2 



Ro(r) 



( dR(r, R ] 
\ dR 



dr z 



Ro(r) 



R(t, R (r)) 2 dn, 



(A6) 



Since, by assumption, R(t = 0) = Rq, 



T = Q 



= 1, this metric is indeed in isotropic form at r = 0. 
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